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We present a new methodology to analyze complicated multi-physics simulations by intro- 
ducing a fictitious parameter. Using the method, we study quantum mechanical aspects of 
an organic molecule in water. The simulation is variationally constructed from the ab initio 
molecular orbital method and the classical statistical mechanics with the fictitious parameter 
representing the coupling strength between solute and solvent. We obtain a number of one- 
electron orbital energies of the solute molecule derived from the Hartree-Fock approximation, 
and eigenvalue-statistical analysis developed in the study of nonintegrable systems is applied to 
them. Based on the results, we analyze localization properties of the electronic wavefunctions 
under the influence of the solvent. 

KEYWORDS: eigenvalue statistics, avoided crossings, quantum chaos, multi-physics simulation, variational 
construction of coupled simulation, electronic states of bio-molecules, one-electron orbital 
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The multi-physics simulation is one of the powerful 
methods to construct complex simulations representing 
realistic systems. Such a calculation is often constructed 
by combining multiple theories with different scales of de- 
scription based on different approximations, e.g., climate 
simulations by fluid dynamics of the atmosphere sur- 
rounded by various external heat and humidity sources,^ 
quantum materials bound to classical large degrees of 
freedom,^' etc. Since reality and accuracy are required, 
those simulations have been larger and more complicated 
year by year. Then, guaranteeing their convergence and 
reliability becomes a more difficult problem. 

In this letter, we study electronic states of a peptide 
in water obtained by a multi-physics simulation which 
consists of ab initio molecular orbital (MO) methods 
and classical statistical mechanics. The quantum chem- 
ical nature of large molecules in a solvent is also one 
of the most popular topics in physics and chemistry. 
In the present calculation, the simulation is constructed 
from a combination of the self-consistent field (SCF) cal- 
culation under the Hartree-Fock (HF) approximation"*^ 
and the statistical mechanics calculation by three di- 
mensional version of the reference interaction site model 
(3D-RISM).^ Although the coupled simulation of RISM 
and SCF itself was developed in 1990's® and has been 
applied to several molecular systems, it attracts pub- 
lic attention again due to the recent developments of 
the distributed computing environments.'' Since the 3D- 
RISM/SCF coupled simulation^ heterogeneously couples 
the different theories both of which require large compu- 
tational resources, the execution of the simulation for 
large solute molecules is difficult not only in an arrange- 
ment of the computer resources but also in theoretical 
aspects of computation such as stability of the conver- 
gence and reliability of the simulation. 

The aim of this letter is two fold: one is to estab- 
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lish the methodology to ensure a proper execution of 
the complicated simulation, and the other is to intro- 
duce the new analysis for a electronic state of molecules 
using eigenvalue statistics.^ The target of our analysis is 
one-electron orbital energies Ej and its wavefunction 
obtained as spatially independent wavefunctions under 
HF approximation. Although Sj and \4>j) are introduced 
through the variational minimization of the total energy, 
they still retain physical reality without Koopman's the- 
orem if we consider energies and wavefunctions of quasi- 
particles and quasiholes.^° Our key to accomplish these 
aims is an introduction of a fictitious parameter which 
can be used not only to investigate the stability and con- 
tinuity of the results but also to perform sophisticated 
analyses^^""*^ based on the eigenvalue statistics. 

At first, we construct 3D-RISM/SCF coupled simu- 
lation according to the standard variational way used 
in ID-RISM/MCSCF." In the present calculation, the 
electronic states of the peptide are obtained by SCF"* 
under the restricted HF approximation (RHF), i.e., all 
electrons are assumed to configure closed orbitals. The 
variational functional is chosen as a total Helmholtz free 
energy ' 

A{X) = £;soiutc({|<^,)}) + Am({|0,)}; A), (1) 

where i?soiuto is an energy of the solute molecule and A/i 
is the excess chemical potential of the solute obtained in 
3D-RISM with the coupling strength A. Statistical cor- 
relation functions of the solvent are calculated under the 
charge distribution of the solute represented by a set of 
the orbitals {|(/>j)}. The variation of the first term in r.h.s 
of (I) with respect to \(j)j) gives usual SCF equations for 
the MO theory while the solvation effect is introduced 
through the second term. Variation with constraints us- 
ing Lagrange multipliers gives conditions to be satisfied. 
Thus, the total free energy (I) of the coupled system is 
stabiHzed by solving these conditions. 

One of the merits of the variational construction of 
complicated simulations is that a set of equations to be 



2 J. Phys. Soc. Jpn. 



Letter 



Author Name 



solved is given by a simple calculus once we define a func- 
tional of the whole system. In the present case, we only 
have to define the solute-solvent interaction potential 



Uj{r; A) = A 



9a ^7 



+ 4e„7<!l — 1 -I — 



(2) 

where is a point charge on a solute atom a, is a 
partial charge on a solvent site 7, ra = |Ra — r| is a 
distance between the solute atom a and the solvent site 
7, Ca-y and 0-07 is given by the standard parameter sets 
of Lennard- Jones potentials. We adopt the method of 
Mulliken population among several schemes to define 
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where C^j are coefficients in a basis-set expansion of 

\<t>3) = E \^-)^-^ (4) 

by the atomic basis function \xv)- This representation re- 
duces computational costs of the whole simulation, and 
can be used even in large-scale simulations for proteins. 
Although the reduction of the cost is often a result of 
compensation in accuracy, it can be retained by combin- 
ing other schemes of the potential representation.^^ 

After the variation with respect to C^j, the solvated 
Fock matrix element F, 

where the overlap integral S^i, is defined by {Xn\Xv)- Vixv 
represents the environmental potential term induced by 
the partial charge distribution of water molecules, 



between /i and v is given by 
= - AF^,5^,, (5) 
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Q(r;A)rfr, (6) 



where R^^ and R^ represent atomic centers of the basis 
functions |x^) and \Xv)-i respectively, and 



A) = ^p^q-yg^{v\ A) 



(7) 



is the partial charge distribution calculated by summing 
partial charges given by the distribution function (77 (r; A) 
over all the solvent sites 7. We can numerically solve the 
Roothaan equation^ 



F(A)C = SCe 



(8) 



to obtain coefficients C^^ of the basis-set expansion as 
elements of the matrix C and the orbital energies Sj as 
diagonal elements of e. 

The actual procedure to execute the coupled simula- 
tion is the following. The SCF calculation is performed 
from an appropriate initial distribution of the solvent. 
For the given distribution of partial charges of the sol- 
vent, the electron density is obtained after a convergence 
of the SCF. Then, 3D-RISM calculation to obtain 5(7 (r) 
is carried out using the charge distribution of the solute 
molecule. Converged 57 (r) determines the environmental 
term W^i, in the solvated-Fock matrix by which the SCF 
calculation is performed again. This procedure repeats 
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Fig. 1. The integrated density of one-electron orbitals of a met- 
enkephalin. Two curves (vacuum and solvated) axe plotted 
around the area near the orbitals of HOMO and LUMO. 



until the whole distribution is converged. From the con- 
verged simulation, we obtain the solute electronic states 
and the solvent distribution functions. 

Before entering the analysis using A, we show results 
of a methionine-enkephalin in two specific cases. This 
peptide consists of a sequence of 5 residues, Tyrosine, 
Glycine, Glycine, Phenylalanine, and Methionine, with 
75 atoms in total. From the standard SCF calculation, 
which corresponds to A = 0, the orbital energies of the 
methionine-enkephalin arc obtained in vacuum, while the 
fuUy-solvated case A = 1 gives the orbital energies in 
aqueous solution. For the electronic ground state, we as- 
sume that 304 electrons occupy 152 independent spatial 
orbitals under RHF. Among several band- like structures 
in the orbital energies, we concentrate the analysis of 
the band just below the highest occupied molecular or- 
bital (HOMO) and the lowest unoccupied molecular or- 
bital (LUMO), which contains 108 orbitals from 45th to 
152th ones. Orbitals strongly localized on Is atomic or- 
bital of C, N, O, and Is, 2s, 2p orbitals of S (known as 
core orbitals) can be found in the region far from the 
band we analyze. Since they exhibit almost no interac- 
tion with other orbitals, we do not consider these orbitals 
in the present analyses. The integrated density of states, 
-O(^) = — Ej), is shown for the band around 

HOMO and LUMO in Fig. 1. In this case, the 152nd 
orbital is HOMO and the 153rd is LUMO. 

Despite that these orbital energies calculated from the 
generalized eigen-equation (8) are not usual eigenvalues 
defined in an autonomous system, they have almost the 
same property as the usual eigenvalues and eigenstates, 
i.e., the orbital energies are real, and the orbital states 
are orthogonal to the state with a different orbital en- 
ergy. Moreover, repulsive interactions between them have 
the same origin as those of the eigenvalues.^" Thus, we 
can perform the eigenvalue-statistical analysis of these 
values after a certain unfolding procedure. In Figs. 2 (a) 
and (b), we show the integrated nearest-neighbor spacing 
distribution obtained after the standard procedure of un- 
folding the spectrum.^ It is known that the fully chaotic 
system shows the Wigner distribution while the regular 
system shows the Poisson distribution. For an interme- 
diate case, we can compare the distribution to a Brody 
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Fig. 2. The integrated distribution of the nearest-neighbor spacing of the orbital energies of a met-enkephalin (a) in vacuum and (b) 
in water. We also show histograms for the nearest-neighbor spacing in the insets, where curves show the Brody distribution (see text) 
with (a) q = 0.52 and (b) q = 0.33. 



distribution. In figs. 2 (a) and (b), we also plot curves for 
the Brody distribution which is obtained by fitting the 
distribution to the integrated Brody distribution, 



PB(s) = l-exp(-as«+^) 
where q represents the Brody parameter, and 



a = 



1 



19+1 



(9) 



(10) 



is obtained by normalization condition. The values of the 
Brody parameters obtained are 5 = 0.52 in vacuum and 
q = 0.33 in water. The smaller value of q reflects the 
smaller interaction between orbitals. 

In order to study further quantum aspects of the solva- 
tion effect, wc enter the detailed analysis of energies and 
wavefunctions to the fictitious parameter A. For chaotic 
quantum systems, the Brody parameter q is usually re- 
lated to the area of stochastic region in phase space. 
However, we do not intend to study any stochastic dy- 
namics behind the electronic states. Instead, we concen- 
trate on the interaction between the orbitals when we 
slightly vary the parameter A, i.e., avoided crossings^'^ 
between the orbital energies. 

The Fock matrix (5) depends on the fictitious parame- 
ter A nonlinearly since A is also included in V^i,, which is 
different from linearly perturbed systems. Even then, we 
can plot the orbital energies to the parameter as usual 
if we perform SCF calculations from A = (vacuum) 
to A = 1 (fully solvated). Figure 3 is a plot of the or- 
bital energies with respect to the parameter A. where we 
chose A^ instead of A as the abscissa because of the non- 
linear dependence on A. Various sizes of avoided crossing 
arc found in the figure, and there seems to be a small 
number of special orbitals which have almost no interac- 
tion with others. Such a separation of orbital states into 
groups without interaction each other can lead to smaller 
values of the Brody parameter. 

The separation of eigenfunctions is often seen in the 
intermediate sciniiclassical limit of quantum chaotic sys- 
tems, i.e., statistical properties in a range around a fi- 
nite energy. A stadium billiard system has a series of 
eigcnfimctions corresponding to a bouncing-ball orbit by 
which the nearest-neighbor spacing deviates from the 



ideal case of random matrix systems. The eigenstates 
corresponding to the bouncing-ball orbit is extremely lo- 
calized in momentum space. ^"'^'^^ In the present results, 
such localization can be observed in the coordinate space. 
Actually, the core orbitals, which we have removed from 
the Brody analysis, are strongly localized on a certain 
atomic basis function. Although the orbitals in the region 
we study here do not localize so strong, it is reasonable 
to analyze the property of orbital states with respect to 
the localization on the atomic basis functions. 

In order to analyze the extent of localization, we define 
the degree of delocalization of j-th orbital state by 



Wi 



(11) 



through the basis-set expansion coefficient C^j in (4). 
This is a representation-dependent quantity, and we de- 
fine Wj on the standard atom-centered basis set {|xi/)} 
in the present work. Although {|Xi/)} is not orthogonal 
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Fig. 3. One-electron orbital energies near the highest occupied 
molecular orbital. The abscissa is A^, the square of the coupling 
strength between the solute and solvent. 
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Fig. 4. Distribution P{Wj) of the delocalization index Wj from 
45th to 152nd orbital states. 



each other and we cannot relate the coefficients C^j di- 
rectly to the extent of the spatial distribution, the value 
of (11), the approximate number of basis fimctions con- 
tained, still represents the relative strength of the delo- 
calization. In Fig. 4, we show the distribution of (11). It 
can be seen that the distribution has two peaks over the 
whole values of A while the places of the peaks slightly 
change. That is, the orbitals are divided into two groups: 
relatively localized orbitals (0 < Wj < 60) and delocal- 
ized orbitals (60 < Wj < 130). This separation of orbitals 
is considered to be the reason why the Brody parameter 
exhibits the intermediate values. 

Before summarizing this work, we give some comments 
for future studies. Since our analysis is based on the 
parameter variation,^"' more sophisticated analysis such 
as the level-curvature analysis^^' and distribution 
of gap sizes of avoided crossings^^' can be performed. 
Since these properties are known to be sensitive for in- 
termediate dynamical systems, it will be more appropri- 
ate for analyzing molecular systems. By the use of these 
analysis, it is desirable to investigate quantum aspects 
of large molecules in water. When we investigate large 
systems such as proteins or nucleic acids, more effective 
methods must be; necidcd. One of such methods is the 
fragment MO calculation^^ by which the orbital energies 
can be obtained with small errors. 

Another interesting problem is whether the regularity 
of the quantum many-body system is influenced by sol- 
vation effects or not, while in the present work we have 
concentrated on the localization of molecular orbitals. 
Although we could not give a clear answer within our 
Brody analysis, such analysis on the regularity will be 
performed through the detailed analysis of nodal pat- 
terns, Husimi representation, etc. of the orbital wave- 
functions with respect to the fictitious parameter. Intro- 
ducing the parameter enables us to assign each solvated 
orbital in terms of the orbitals in vacuum. This proce- 
dure will successfully be done by connecting characteris- 
tic features of the orbital over avoided crossings since it is 
possible even in strongly chaotic systems^^ if many pair- 
wise avoided crossings are found. Prom a technical point 
of view, the continuous parameter is also important as a 
computational technique. In our case, the stability of the 
result is numerically supported by varying the parame- 
ter slightly, since physically meaningful quantities can 



only be obtained as stable results under small pertur- 
bations. In order to obtain stable and significant results 
from complicated realistic simulations, the fictitious pa- 
rameter approach will be one of the powerful methods. 

In summary, we have studied statistics of the orbital 
energies and localization of orbital states through the 
variation of the fictitious parameter. This is our first re- 
port on the eigenvalue-statistical analysis of the orbitals 
as a new approach to quantum many-body systems. We 
hope that these new methods will be used to study quan- 
tum aspects of large bio-molecules as well as to complete 
large-scale realistic simulations. 
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